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Abstract. Nonlinear effects are crucial in order to compute the cosmological matter 
power spectrum to the accuracy required by future generation surveys. Here, a 
new approach is presented, in which the power spectrum, the bispectrum and higher 
order correlations, are obtained - at any redshift and for any momentum scale - by 
integrating a system of differential equations. The method is similar to the familiar 
BBGKY hierarchy. Truncating at the level of the trispectrum, the solution of the 
equations corresponds to the summation of an infinite class of perturbative corrections. 
Compared to other rcsummation frameworks, the scheme discussed here is particularly 
suited to cosmologies other than ACDM, such as those based on modifications of 
gravity and those containing massive neutrinos. As a first application, we compute 
the Baryonic Acoustic Oscillation feature of the power spectrum, and compare the 
results with perturbation theory, the halo model, and N-body simulations. The density- 
velocity and velocity-velocity power spectra are also computed, showing that they are 
much less contaminated by nonlinearities than the density-density one. The approach 
can be seen as a particular formulation of the rcnormalization group, in which time is 
the flow parameter. 
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1. Introduction 

Future generation galaxy surveys are going to measure the statistical properties 
of matter distribution to an unprecedented accuracy, providing informations on 
fundamental questions such as the nature of Dark Energy (DE) [U Ej and the absolute 
scale of neutrino masses (see, for instance, [3], and references therein). In particular, the 
location and amplitude of Baryon Acoustic Oscillations (BAO), wiggles in the matter 
power-spectrum produced by the coupling of baryons to radiation by Thomson scattering 
in the early universe, for wavenumbers in the range k ~ 0.05 - 0.25 /iMpc" 1 , have the 
potential to constrain the expansion history of the Universe and the nature of the DE. 
BAO's have recently been detected both in the 2dF and SDSS surveys data [U (5J El [7] , 
and are going to be measured in the near future in a series of high-redshift surveys (see, 
for instance [HIE]). 

A reliable comparison between theoretical models and observations requires going 
beyond the linear order in perturbation theory [101 [HI [EJ [T3], [TH [15]. The more 
established way to deal with nonlinearities is by means of A^-body simulations (for recent 
applications to the matter power spectrum (PS) see for instance [TBI [TT1 ITS] ). However, 
in order to gain the required sensitivity, very large volumes and high resolutions are 
required, with the consequence that, due to time limitations, only the 'vanilla' type 
ACDM cosmologies have been investigated so far. Fitting functions for the nonlinear 
PS, such as those based on the Halo model e.g. refs. [191 120]) are uncapable to reach 
the required level of accuracy [TB] . 

An alternative, semi-analytical approach is represented by perturbation theory 
(PT) (for a review, see [21]), which has recently experienced a renewed interest, mainly 
motivated by two reasons. First, next generation galaxy surveys are going to measure 
the PS at large redshift, where the fluctuations are still in the linear regime and 1-loop 
PT is expected to work [TQl [Tl] . Second, techniques for the resummation of some classes 
of perturbative corrections to all orders have been recently developed. These techniques, 
based on field theory tools such as Feynman diagrams and the renormalization group 
(RG) [22j [231 E2 EH] have been shown to be able to render a PS in agreement with 
that from A^-body simulations down to z = [131 HH US]- F° r related work, see 

[SlESlEBlETlEHlEniEolEI]. 

While these semi-analytic models have proved themselves a viable alternative to 
A^-body simulations (at least for what concerns the Dark Matter (DM) PS), they still 
suffer from some limitations. First of all, these methods are formulated for an Einstein- 
deSitter (EdS, Q m = 1) cosmology, and further extended to more general ones (e.g. 
ACDM) by replacing the EdS linear growth factor for the growing mode, Df dS = a, with 
that of the other cosmology, D + (a), where a is the scale factor. While this procedure 
is correct at the linear order, it introduces some inaccuracies at higher orders as soon 
as the condition fi m //J = 1 (with / + = d In D + /d In a) fails [32]. The problem is 
that in schemes such as [22l [23l HSl [14] there is no way to independently assess the 
validity of such approximation. Moreover, and more importantly, such approaches 
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cannot be trivially extended to cosmologies in which the linear growth factor is also 
scale dependent, D + = D + (k, a), as is the case, for instance, when massive neutrinos 
[3] contribute to the DM. 

On a more technical ground, while the leading corrections for the PS where correctly 
identified and computed in [131 EH [IS], for the bispectrum (BS) the task turns out to 
be much more involved (for a computation in a different approach, see [33J ) . Finally, a 
systematic scheme of approximations was not univocally identified in such frameworks, 
and the independent assessment of the size of the corrections one is neglecting is not a 
straightforward task. 

In this paper we present an alternative approach to the resummation of infinite 
classes of PT corrections. The method is based on a straightforward application of 
the 'equations of motion', i.e. the continuity, Euler, and Poisson equations, to the 
computation of correlation functions (PS, BS, etc.) for the density and the velocity 
fields. The correlation functions evolve in time according to a (truncated) system 
of differential equations, quite similar to the familiar BBGKY hierarchy ones [34]. 
Compared to the approaches discussed in [221 [231 [131 EH], the one we are going to 
present is much more free from field theory technicalities and, as such, will be more 
easily implemented by people without a specific field theoretical background. The 
correlation functions obey a system of integro-differential equations, closed by a well 
defined approximation procedure, whose solution gives the PS, BS, etc. at different 
times and for any momentum configuration. 

Remarkably, as we will discuss, the present formulation is much more flexible than 
the previous ones, and can be straightforwardly applied to more general cosmologies, as 
those based on modifications of gravity and those containing massive neutrinos [35] . 

As a first illustration of the possible applications of the method, we compute 
the PS in the BAO region. We discuss the simplest non-trivial (i.e. beyond 1-loop) 
approximation level and show which class of PT corrections it corresponds to. We 
consider two cosmologies, a 'vanilla' ACDM, and one in which DE is a quintessence 
with equation of state w = —0.8. For the ACDM case, we compare our results with 
those obtained by ./V-body simulations [10] and by other approaches, and for both cases 
we assess the error made by not properly treating the decaying mode. We find that it 
can be as large as ~ 1% only for z close to zero and scales smaller than those relevant 
for the BAO. 

The density- velocity and velocity- velocity power spectra are also computed, showing 
that they are much less contaminated by nonlinearities than the density-density one. 

The paper is organized as follows. In Section [2] we discuss the fluid equations for 
a large class of cosmologies and derive, starting from those, the evolution equations 
for the PS's and the BS's. We also state the only approximation made in this paper, 
consisting in neglecting the effect of the evolution of the connected four point functions, 
the trispectrum. In Section [3] we derive a formal solution of the system of equations 
in our approximation and show to which class of PT corrections it corresponds to. In 
Section @] we compare the present approach to those of refs. [221 [231 EH [H], and also 
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with the one, closer in spirit, of ref. In Section Owe give and discuss our numerical 
results and compare it to other approaches. Finally, in Section [6] we summarize our 
results and discuss future lines of work. 



2. Nonlinear fluid equations in general cosmologies 

Our starting point is given by the fluid equations, 

^ + v-[(i + <Uv] = o, 

or 

-^ + H(v+[Av]) +(v-V)v = -V0, 
or 

V 2 ( p = ^H 2 Q m (6 m +[B5 m }), (1) 

which are the continuity, Euler, and Poisson equations, respectively, generalized, with 
respect to the EdS case, by the presence of the functions fl m (r), *4(x, r), and £>(x, r). 
We will work in conformal time, r, and will denote by 7i = dloga/dr the Hubble 
parameter, while <5 m (x, r) and v(x, r) will be the DM number-density fluctuation and 
the DM peculiar velocity field, respectively. The square brackets indicates a convolution 
in real space, i.e., 

[0 x] (x, r) J -0^ 0(x - y , r) X (y, r) . (2) 

The functions A(x, r), £>(x, r), and Q m (r), parametrize different cosmologies. In 
the following, we list some interesting cases; 

1) A = 0, B = 0, and fi m (r) = 1, corresponds to Einstein-de Sitter cosmology; 

2) A — 0, B — 0, and Vt m {r) = 1 + pg/ ' p° m f{r) , corresponds to a flat cosmology 
containing DM and some non-clusterized quintessence fluid. The case of constant 
equation of state is given by /(r) = a(r)~ 3w , and ACDM corresponds to w = — 1; 

3) Non- vanishing values for both A and B occur when particles' geodesies are modified, 
as it happens for instance when the DM particles couple to some light scalar 
field [361 EH EE] or in Scalar- Tensor theories of gravity [391 HO]- I n the latter 
case, if a(ip) represents the coupling between matter and the scalar field (in 
terms of the Brans-Dicke parameter u, we have a 2 = l/(2u + 3)) one has 
*4.(k, r) = A(t) = ad(p/d\og a and £>(k, r) = 2a 2 (l +m 2 J a 2 /k 2 ^ , where 

is the scalar field mass, and we have transformed to Fourier space. 

The approach we are going to present in this paper can be straightforwardly extended to 
the case of many fluctuating components. In this case, the RHS of the Poisson equation 
becomes 

!N 2 E^ (3) 

A i 

where the sum extends to all the components having non- vanishing density fluctuations. 
For each component, nonlinear continuity and Euler equations should be taken into 
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«(q,p) = - 2 — , /%,p) = — 7T-2~2 — ' ( 6 ) 



account. In some cases, however, the non-DM components have fluctuations well inside 
the linear regime in all scales of interest. This happens, for instance, for neutrinos and 
for most quintessence models. In these cases, one can approximate eq. ([3]) as 

where the terms inside the parenthesis are taken from linear perturbation theory. Then, 
in this approximation, the system reduces to eqs. ([T|) for a single DM component with 
A = and B(k, t) = Ei ft< Sf n /(n m 

Defining, as usual, the velocity divergence #(x, r) = V • v(x, r), and going to 
Fourier space, the equations in (JT]) give 

or 

+ J d 3 qd 3 p5 D (k- q-p)a(q,p)9(q,T)5 rn (p,r) = 0, 

d9 ^ T) +H(1 + A(k, T))9(k, r) + ^H 2 (1 + B(k, r))Q m (r)S m (k, t) 

+ |rf 3 qd 3 pfe(k-q-p)/5(q,p)^(q,r)0(p,r) = O. (5) 
The nonlinearity and non-locality of the Vlasov equation survive in the two functions 

(p + q)-q „/ v (p + q) 2 p-q 

which couple different modes of density and velocity fluctuations. Notice that the 
functions A and B only appear in linear terms. 

One can write Eqs. (jSJ) in a compact form. First, we introduce the doublet <p a 
(a = 1,2), given by 

(<Pi(Kv)\ = -r,( tm{k,r}) \ 

\^0cr,)J- { -9(k, V )/H) ' {) 

where the time variable has been replaced by the logarithm of the scale factor, 

i a 

V = log — , 

"'in 

ai n being the scale factor at a conveniently remote epoch, were all the relevant scales 
are well inside the linear regime. 

Then, we define a vertex function, 7 afec (k, p, q) (a, b, c, = 1,2) whose only 
independent, no n- vanishing, elements are 

7m (k, p, q) = - 8 D (k + p + q) a(p, q) , 

7222 (k, p, q) = 5 D (k + p + q)/5(p,q), (8) 

and 7m (k, p, q) = 7m(k, q, p). 

The two equations ([5]) can now be rewritten as 

dr,<p a (k,rj) = -n ab (k, r])(p b (k,r]) + e v -f abc (k, -p, -q}(p b (p,-n)tp c (q,7]),(9) 
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where 



fi(k, V) 



1 -1 



) 



(10) 



Repeated indices are summed over, and integration over momenta q and p is understood. 
Notice that all the effect of changing the cosmology from EdS is contained in the matrix 
Q ab , the vertices being universal. 

The ?7-evolution of the correlation functions can be derived by iterating the 
application of eq. (jHJ), as follows: 



+ e^acd^cipdifib) + e v -f bcd (ip a ip c ip d ) , 

d T , (VaWPc) = - QadiVd^b^c) - ^bdi^a^d^c) ~ ^cd(^a^b^d) 
+ e r 'lade(<Pd<Pe<Pb<Pc) + e^lbdeiVaVdfe'Pc) 

+ e n j cde ((p a (p b (p d (p e ) , 

9 V {fafbfcfd) = ■■ ■ 



In order to have more compact equations, we have omitted the momentum and in- 
dependence of the correlation functions. All the fields are evaluated at the same 77- value. 
Next we express the two -three -and four-point correlators as 



{<p a {k, rj)ip b (q, rj)) = S D (k + q)P a fe(k , rj) , 

(<Pa(k, ??)<Mq, v)<Pc(p, rj)) = M k + q + p)B abc (k, q, p; 77) , 

(<p a (k, rj)ip b (q, rj)ip c (p, rj)ip d (r, rj)) = 

[5 D (k + q) 5 D (p + r)P ab (k , v)Pcd(P , V) 
+ 5 D (k + p) S D (q + r)P ac (k , rj)P bd (q , rj) 
+ 5 D (k + r) 5 D (q + p)P ad (k , r])P bc (q , 77) 



where P ab (k,r]) is the PS, B abc (k, q, p; 77) the BS, and Q a bcd(k , q , p , r , 77), the 
connected part of the four-point function, is the trispectrum. 

In this paper, we will make the approximation Q a bcd = 0. It should be emphasized 
that this choice, although allowing us to split the four-point functions in terms of two- 
point ones as in the Wick theorem, by no means amounts to consider the fields tp a to 
be Gaussian, since the BS is fully taken into account. In the following Section, we will 
show explicitly by using diagrammatic methods, which class of non-linear interactions 
are kept into account by this approximation. 



(11) 



+ 5 D (k + p + q + r) Qabcdik , q , p , r , 77)] 



(12) 
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The first two equations in (ITTj) form now a closed system, given by: 
d v P ab (k , rj) = -fi ac (k , 7])P cb {k , 77) - tt bc (k , ?7)P ac (k , 77) 

+ e v J al 3 q [7 acrf (k, -q, q - k) B bcd (k, -q, q - k; 77) 

+B acd (k, -q, q - k; rj) 7 fecrf (k, -q, q - k)] , 

d v B abc (k, -q, q - k; 77) = -Q a d(k , r])B dbc (k, -q, q - k; 77) 

- ft M (-q , r)) B adc(k, -q, q - k; 77) 

- fi cd (q - k , r])B abd (k, -q, q - k; 77) 

+ [7ad e (k, -q, q - k)P db (q , 77)P ec (k - q , 77) 
+ 7&de(-q, q-k, k)P dc (k - q , 77)P ea (k , 77) 
+ 7ccfe(q - k, k, -q)P da (k , r/)P e6 (q , 77)] . (13) 

Notice that we have restored the momentum dependences, and that the only momentum 
integration is the one explicitly indicated in the first of eqs. (1131) . In Section [5] we 
will present the result of the numerical integration of the equations above, without 
any further approximation, so, the reader interested in the performances of the 
present approach can skip directly to that section. In the next two sections, we 
discuss the relation of this method with PT, with the resummation schemes of refs 
[221 E31 EHl CEH GS], and with the RG approach of ref. [26]. 

3. Analytic Solutions 

The formal solution of the system (fTBl is given by 

P afe (k , 77) = g ac (k , 77, 0) g bd (k , 77, 0)P cd (k , 77 = 0) 

+ J dr]'e v ' J d 3 qg ae (k,r],r]')g bf (k,r],r]') 

x [y ecd (k, -q, q-k) B fcd (k, -q, q-k; 77') 

+ 7/cd(k, -q, q-k) B ecd (k, -q, q-k; 77')] , 



x [idghik, -q, q - k)P e9 (q , rj')P fh (q - k , 7/) 
+ 1e 9 h(-<l, q-k, k)P /ff (q - k , r]')P dh (k , 7/) 
+ 7/5 ,(q-k, k, -q^k^OP^q,//)] , (14) 

where g a b(k , 77, 77') is the linear propagator [221 E31 [131 E] which gives the 77-evolution 
of the field at the linear level, ^(k, 77) = g ab {k , 77, rj')Lp^{k, 77'), where the subscript "L" 



B abc (k, -q, q - k; 77) = 

g ad (k , 77, 0)g be (-q , rj, 0)g cf (q. - k , 77, 0)P de/ (k, -q, q - k; 77 = 0) 
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Figure 1. The 0(7°) contributions to the power spectrum and the bispectrum. 



indicates the linear approximation. In order to check the validity of the solutions above, 
the two following properties of the propagator should be used, 

•9^6 (k, r/,?]') = -n ac (k,r])g cb (k, 77,77') , g ab (k,r),r]) = 5 ab . (15) 

The explicit form of the linear propagator for the general class of cosmologies 



encompassed by eq. ffTUl) is given in Appendix A 
It is instructive to expand eqs 



1141) in powers of the interaction vertex •jabc, m 
order to understand what class of peturbative corrections are mantained, and what are 
neglected, by doing the approximation Q a bcd = 0. The lowest order, corresponding to 
linear perturbation theory, is obtained by setting ^ abc = on the RHS of eqs. ((Hj), 
yelding, 

P^ik , 77) = g ac (k , 77, 0) g bd (k , 77, 0)P cd (k , V = 0) , 



B a bc (k, -q, q - k; rj) = 

9ad(k , 77, 0)g be (-q , 77, 0)# c/ (q - k , 77, 0)B def (k, -q, q - k; 77 = 0) , (16) 

where we explicitly see, besides the initial PS, the initial BS evaluated at 77 = 0. Notice 
that, due to the e _T? factor in front of eq. (ED), the linear PS is independent on 77, on the 
growing mode. 

In the language of Feynman diagrams introduced in [T3j, the linear order solution 
can be represented as in fig. [lj where the lines correspond to propagators, the empty 
box to the primordial PS and the empty triangle to the primordial BS. The next order 
(0(7)) correction for the BS is obtained by inserting P^ h in place of P ab at the RHS 
of the second of eqs. fTl4|) . The result is represented at the bottom of fig. [2j where the 
vertex represent the interaction ^ abc . Inserting the BS at this order (that is, the sum of 
the 0(7°) and 0(7) contributions of figs. [1] and El respectively) in the first of eqs. (fT4j) . 
we get the 0(7 2 ) and 0(7) contributions to the PS, represented in the first two lines of 
fig. [21 At this order, the result for the PS coincides with the standard 1-loop expression 
in eulerian perturbation theory [21], as can be checked by using the composition rule 
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Figure 2. The 0(7) contributions to the bispectrum and the 0(^ 2 ) (first line) and 
0(7) (second line) ones to power spectrum. 



for the propagators, 



g ab (k , 77, r]')g bc (k , 7/, rj") = g ac (k , 77, 77") . 



(17) 



Iterating the procedure one more time, the first differences with respect to perturbation 
theory arise. Considering the BS, 0(7 3 ) and 0(7 2 ) corrections are generated by 
substituting one of the two power spectra lines in fig. [2] with the 0(7 2 ) or 0(7) 
contributions to the PS, one possibility being represented on the left of fig. El However, in 
perturbation theory, there are other 0(7 3 ) corrections which cannot be obtained in this 
way, such as that on the right of fig. In a field theoretic language, as a consequence of 
the approximation Q a bcd = 0, we are neglecting vertex renormalization, while including 
the renormalization of the PS. This finding generalizes at any order. The n-th iteration 
leads to contributions which are given by sums of terms having the same structure as 
fig. [21 with the legs containing a box replaced 0(n — 1) contributions to the PS, while 
the vertices, propagators, and triangles are left untouched. This is well summarized 
by the expressions in eq. (fl4|) . where the PS is formally 1-loop (there is a momentum 
integration, d 3 q), while the BS is formally tree- level (no momentum integration). 

The results presented in Sectj5] correspond to the resummation of this class of 
contributions to all order in 7. 

4. Comparison to other approaches 

From the discussion of the previous section, one can realize that the results of linear 
perturbation theory for the PS are recovered by setting 7 a b c to zero on the RHS of the 
integro-differential equations in (|T3|) . The 1-loop approximation corresponds instead to 
setting P a b(p, v) = Pab(p,V = 0) 011 the RHS of the second of eqs. ( jl3l) . while keeping 
the first one in its full form. 
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Figure 3. 0(j 3 ) corrections to the BS included (left) and not included (right) in the 
approximation scheme considered in this paper, namely Q a bcd = 0. 

Finally, keeping the ^-dependence of the PS at the RHS of both equations 
corresponds to resumming, to all orders, the perturbative corrections in which the 
interaction vertex is kept at its tree-level form, 7 a b c , i.e. in which no diagrams like the 
one at the right of fig. [3] are included. This class of corrections is the same which was 
considered in refs. [Hj, where RG equations were also solved under the approximation 
of keeping the tree-level expression for the vertex. Indeed, the results of the present 
paper compare very well with those obtained via the RG. However, there are also some 
practical differences between the two approaches, which are worth while mentioning. 

First of all, in the RG CcLSG, clS well as in the renormalized perturbation theory of 
Crocce and Scoccimarro [22], [23], the computation of the PS requires two steps, each 
one implying further approximations besides the one of keeping the vertex at tree-level. 
In the first step, the propagator is computed. As was shown in [22} [231 E3 EH] a simple 
analytic expression, exact in the limit of large external momentum, can be obtained 
for the propagator at any order in perturbation theory. However, the use of this large 
momentum limit turns out to give inaccurate results [15], so that subleading corrections 
have to be taken into account. In the second step, the 'resummed' propagator is used to 
compute the PS. The calculation is performed either by making extra approximations 
to the relevant RG equations [14], or at finite loop order, as in |15j . 

The remarkable feature of the present approach is that there is no need to compute 
the propagator in order to get the PS (or the BS, etc.). The information on the time 
evolution, which, in the approaches of refs. [22} [23], [T3l HU [15] is contained in the 
propagator, here is treated exactly by the differential (in time) structure of the equations. 
Once the Q a bcd — truncation is performed, the system is solved without any further 
approximation. Also, the next level of approximation is clearly identified in keeping 
Qabcd and truncating at the five-point correlator level, and so on. 

Moreover, as we already discussed in the Introduction, using differential equations 
in time has the other precious virtue of allowing a clear treatment of the general class 
of cosmologies described by eqs. [U Strictly speaking, approaches such as [13], HU [15] 
are exact only in the Einstein-deSitter case (Qm = 1), in which the matrix 0(k, 77) 
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appearing in the the equation of motion reads 
I 1 - 1 

n = 33 

V 2 2 

In cosmologies such as ACDM, where Q m < 1 and the linear growing mode for matter, 
D + , is scale- independent, if one makes the following replacements 

^^ log 7w — V' ^(M) -> ——— , (19) 
with / + = dlogD + /dloga, one gets equations of the same form as eq.(jn]) but with 

( 1 -M 

n(r/) = _3fi„ 30^ . (20) 
V VI Vl I 

The phenomenologically interesting cases of ACDM and (non-clusterizing) quintessence 
are then treated in [221 1231 EES HU [15] by using eq. (1T9"]) and approximating eq. (1201) 



with ( fl8l) . It is easy to realize (see Appendix A) that in this approximation the growing 



mode of matter perturbations is treated properly, whereas the decaying mode has the 
wrong time dependence and the wrong ratio between density and velocity perturbations. 
Since the vertex 7 a fe c mixes growing and decay modes, this approximation is expected to 
work well in the linear regime, and to fail when nonlineari^ become importanfl The 
frameworks discussed in [221 1231 [T3l [HI EES] have no way of assessing the validity of the 
approximation, while, in the present approach, one can directly test it by comparing 
the results obtained with the constant and the 77- dependent O matrices. In the next 
section we will discuss such a comparison (see fig. [7]). The most generic case of scale 
dependent growth factors is precluded to the approaches such as [221 [231 [131 CHI [15] - 
while it can be properly treated here. The most relevant example of massive neutrinos 
will be the subject of a future publication [35] . 

Evolution equations for the PS depending on time (or, better, on the growth factor), 
were discussed by Mc Donald in [26]. They are equivalent (at least in the approximated 
approach to ACDM discussed above) to make the approximation 

Pab(p, rj) ~ u a u b P(p, 77) , u =^j^J, (21) 

in the last of eqs. fll4l) . and in inserting the BS obtained in this way into the first of 
eqs. (fTBl . By doing so, one obtains an evolution equation whose RHS has the structure 
of the 1-loop correction to the PS, in which the linear (//-independent) PS has been 
replaced by the renormalized (^-dependent) one. As we will see in the next section 
(see figs. El EI), the different components of the PS exhibit sizably different evolutions, 
therefore one expects that the approximation in eq. ([2"T]) would lead to results far from 
the accuracy required in applications such as the BAO, as it seems to be the case of 
ref. [261. 



\ Another way of seeing the problem [32] is to notice that if fl m / f+ 7^ 1 the equations of motion (J9j) 
do not admit separable solutions of the form </? a (k, i]) = J2 n -D+( 7 ?)Xa(k). 
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Figure 4. Power spectra at redshift z = 1 (divided by a smooth one). The continuous 
line is the result of the present paper, compared with linear theory (dotted), 1-loop 
PT (dash-dotted), the halo approach of ref. ^2U\ (dashed). The dots with error bars 
are taken from the TV-body simulationd of ref. [TO]. The background cosmology is a 
spatially flat ACDM model with n° A = 0.73, = 0.043, h = 0.7, n s = 1, <r 8 = 0.8. 




Figure 5. The density-density (Pss), density- velocity (Pse), and velocity- velocity 
(Pee) power spectra at z = 1. The dotted line is linear theory, and the dash-dotted 
one the linear theory multiplied by exp(— fc 2 (7 2 (e' , — l) 2 ). 



Flowing with Time: a New Approach to Nonlinear Cosmological Perturbations 



13 




Figure 6. Same as fig. El but at z = 0. 



5. Numerical Results 

In this section, we present the results of a numerical integration of the equations ([[3 



Our code is built along the lines described in|Appendix B| As benchmark cosmologies 



we will consider a fiat ACDM model, with tt° A = 0.73, tt° b = 0.043, h = 0.7, n s = 1, 
as = 0.8, and a quintessence model with equation of state w = —0.8 and the other 
parameters set at the same values. The linear growth factors are scale-independent in 
both cases, so we will perform the change of variables in eq. ( Tl9l) and use the exact 
-//-dependent - matrix ft in eq. (1201) . 

The integration starts at a rji n = corresponding to a high redshift. In this paper, 
we set z in = 35. The initial conditions at rj = are the linear PS, as obtained by running 
the CAMB code [UJ, and vanishing BS, i.e., we are ignoring all non-gaussianities 
generated at higher redshifts. Then, the evolution equations are evolved down to low 
redshift, giving the PS and the BS for any momentum scale. In fig]U we give the 
(density-density) PS for our ACDM cosmology, divided by a smooth spectrum, obtained 
by setting = 0. The results of this paper are represented by the continuous line. For 
comparison, we also give the results of linear theory (dotted) of 1-loop PT, dash-dotted, 
of the halo approach of ref. [20] (dashed), and of the A^-body simulations of ref. [TO] . 
Compared to the RG approach of refs. [T3l [14"] . the agreement with A^-body simulations 
is equally good in the BAO region, but here, it extends to higher fc's. 

In figs. [5] and [6] we plot the density- density, density-velocity, and the velocity- 
velocity PS, given by Pn, Pu, and P22, respectively. As it was discussed in [TBI [T4l ITS"] . 
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the exact PS has the structure 



P ab (k, V ) = G ac {k, V , 0)G bd (k, rj, 0)P cd (k, 0) + P%{h, rj) 



(22) 



where G ac (k, rj, 0) is the full propagator. The first term on the RHS can be approximated 
as PI 



the velocity dispersion. This contribution is also plotted in the figures, with a dash- 
dotted line. In this approximation, the first contribution to eq. ( J22l) is the same for the 
three components of the PS. The difference between them is then entirely due to the 
mode-coupling contribution P^{k, rj), which has the effect of smoothing out the BAO 
feature [IB] . As we see from the figures, the velocity- velocity and the density- velocity 
PS are less affected by P n than the density-density one, therefore they should be taken 
into serious consideration as alternative keys to measuring the BAOs. 

Finally, we address the issue of the error made by approximating the exact fi in 
eq. ( 1201) with the constant one in eq. (TP8|) . as is usually done in PT as well as in the 
approaches of refs.[22], [231 Efl EH EE]- This approximation should become less and less 
accurate at decreasing redshifts. Also, it should have some dependence on the cosmology 
under consideration, since the ratio Q m /f+ increases from 1 at high redshift to 1.157 at 
redshift zero for ACDM and to 1.179 for the w = —0.8 quintessence cosmology. In fig. [7] 
we plot the relative difference between the PS computed with the time-dependent matrix 
n and the approximated, constant one. Continuous lines are for ACDM (w = —1), and 
dashed ones for a quintessence with constant equation of state w = —0.8. For each 
cosmology, the lines from bottom to top correspond to z — 1, 0, respectively. We see 
that, in both cases, the relative error in the BAO range of stays below the per mill level 
at z = 1, and below the percent at z = 0. Therefore, for high redshift surveys, our 
results indicate that the approximation is well motivated. 

6. Summary 

Compared to previous semi-analytic approaches to the resummation of perturbative 
corrections, the scheme presented in this paper has some clear advantages. First of 
all, it can be easily extended to a large class of cosmologies, among which the case of 
massive neutrinos and of scalar-tensor modifications of gravity are the most notable 
ones. In practice, the approach in its present form can be used for all models in which 
the nonlinear terms in the continuity and in the Euler equations are the same as for the 
EdS case. The evolution can be seen as a series of time steps (in the limit of vanishing 
step length), during which the fields evolve with the -possibly scale-dependent- linear 
growth factors. At the end of each step, the growing and decaying mode are mixed by 



with P L 



e (- fe ^(e,-l)3)pL (A . Q) 

(k, 0) the linear PS, and 



(23) 




(24) 
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Figure 7. Relative difference between the power spectrum computed with the 
time-dependent matrix f2 in eq. (|20[) and the approximated, constant one, cq. (I18D . 
Continuous lines are for ACDM (w = —1), and dashed ones for a quintessence with 
constant equation of state w = —0.8. For each cosmology, the bottom (top) lines 
correspond to z = 1 (z = 0). 



the universal nonlinear terms. Therefore, the different cosmologies enter only via the 
different background evolution and the different linear growth factor, all the information 
being contained in the f2 matrix of eq. ( TTOT) . 

The second advantage is the possibility of formulating a clean and systematic 
scheme of approximations. The next level, compared to the one considered here, is 
obviously given by the inclusion of the trispectrum in the system of equations, which 
will correspond to including vertex corrections. Judging from the present agreement 
between our results and TV-body simulations in the BAO region, one would infer that 
such improved approximation will mainly affect higher momentum scales. 

For a more practical point of view, though the present approach is completely 
compatible with those of refs. [221 [221 [EH [H] (see the discussion in Section @J, here field 
theory tools such as Feynman diagrams, generating functionals, and RG flows, are not 
necessary in order to compute, for instance, the PS and the BS. What is needed, at 
the end, is just a code to solve the differential equations in eq.l fTBl . Therefore, a field 
theoretic background is not a prerequisite to use this method. 

The results we obtain are in good agreement with the N-body simulations of ref. 
[TO] , up to k of order 0.3 — 0.35 h/Mpc, showing that the method is able to describe (from 
first principles) the non-linearities due to mode-mode coupling in the BAO region. We 
also showed that the velocity-density and the velocity-velocity PS's are less contaminated 
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by the mode-mode coupling and, since they should also be much less contaminated by 
bias, their possible role in the future BAO measurements should be studied in more 
detail. 

Finally, we showed that the role of the decaying mode can be safely ignored for 
z ~ 1, and can lead to percent-level effects only at lower redshifts and scales higher 
than those relevant for the BAO's. 

The importance of nonlinear effects for the determination of the upper bound on the 
neutrino mass scale was emphasized in ref. [32], where 1-loop PT was used to compute 
the PS. A better nonlinear approach is clearly desirable to derive firmer neutrino mass 
bounds [35], possibly including the running of the trispectrum, in order to reach higher 
momentum scales. 



Appendix A. The linear propagator in general cosmologies 

The linear equation 

d v ip a (k,7]) = -fi o6 (k, r))<fi b (k,r)), (A.l) 
has solutions of the form 

{nL)) M (A - 2) 

where (p and / satisfy, 

d v <p= -(fin + Qnf)<P 

= n«(/ -/+)(/-/-), (A.3) 

with 



_ /t fi 22 - fi n T J(tt 22 - fi n ) 2 + 4fi 12 fi 21 

/±(k, V) = " — • ( A ' 4 ) 

Solving the equations flA.31) . we get the evolution in rj of the upper and lower components 

of dS2D, 

f(vMv) = e-^ dxin2l/f+n22) f(v'Mv') 

= e-^ dx{Qll+Ql2f) lplf( V 'M V ') . (A.5) 

fw) 

Next, we identify two independent solutions of eq. (1A.3[) . which will serve as basis for 
the expansion of the most generic solution of eq. flA.lj) . In most cases, it is convenient to 



choose as basis solutions the growing and decaying modes. For instance, if the Universe 
at high redshift tends to Einstein-deSitter, one can identify the basis solutions by setting 
their initial conditions as 

/±(M<) = /±(Mi) , (A.6) 
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at some early time rji, with f± given in eq. (lA.4j) . 

We define the instantaneous projectors on the two basis solutions, for each 77 > rji, 

/+(M) ) = ( ) ' 



M + {k,T]) 



which is given explicitely by 

M+ = T±h ( L -n ) ' (A ' 8) 

The projector on the other mode is given by M~ = 1 — M + . 

The propagator, defined as the operator giving the linear evolution of the field tp a , 
i.e, <p a (k,r]) = g ab (k,r],rf)ip h (k,rf), is then given by 

g(M, rf) = e -^^M 1 1 \ M+M) 

+ e -/>(a 11+ o 12 /_) (1 M - Mh (A . 9) 

\ f-(Kv') J 

for rj > rf , and zero otherwise. One can check that it satisfies the properties 

d v g(k , rj, rf) = - ft(k , 77) ■ g(k , 77, 77') , 
lim g(k, 77,77') = 1 , 

7]/—>T] 

g(k, 77, rf) ■ g(k, 77', V ") = g(k, 77, 77") , (A. 10) 

and that, in the Einstein-deSitter case (flu = —Qu = 1, ^22 = —^21 = |), it reduces 
to the propagator considered in [Hj, with / + = 1, /_ = —3/2. 

Appendix B. Putting the equations in a more numerical-friendly form 

In case of /c-independent (but 77-dependent) Q matrices, the evolution equation for the 
PS, eq. (fTBl can be rewritten as 

d v p ab (k) = -n ac p cb (k) - n bc p ac (k) 

Air f°° [i r_ ~ 

+ e ri — dqq dpp % cd (k,q,p) B bcd (k,q,p) 

K Jk/2 J\a-k\ 1 



-B acd (k,q, p)%cd(k,q,p) , (B.l 



where the 77-dependence of the PS, the BS and the fl matrix is understood, and we have 
defined 

%bc(k, q, p) = 7abc(k, q, p) | p= _ (k+q) , (B.2) 

and analogously for B a b c (k,q,p). In numerical applications, the integrals appearing in 
the above expression are quite time consuming to perform, both because the integration 
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interval for the internal variable (p) depends on the external one (q), and because the 
integrand contain a function, B be f(k, q,p), which has to be computed via an independent 
evolution equation, eq. ( |T3|) . for each set of its three variables. In order to put the 
equations in a more manageable form, it is convenient to define the following quantities, 

Iacd,bef(k) = / dqq / dpp- % c d(k,q,p) B bef (k,q,p) + (q <h> p) . 

Jk/2 J\q—k\ 6 L J 

(B.3) 

Then, eq. (IB.ip can be written as 

d v P a b(k) = -Vt ac P cb {k) - n bc P ac (k) 

47T 

+ e* 1 — [I acdMd (k) + I bcd , acd (k)} . (B.4) 

The evolution of the I a cd,bef(k) can be deduced from that for the bispectra, i.e. from 
the second of eqs. (fT3l . and is given by 

9r] Iacd,bef(k) = —^bgl acd,gef(k) — Q eg l a cd,bg f(k) 

- Qfgl a cd,beg(k) + 2e v A acdMf (k) , (B.5) 
with the integrals A acdMf (k) given by 

f°° f q 1 
A ac d,bef(k) = dqq dpp - {% cd (k, q,p) [% h (k, q,p)P ge (q)P hf (p) 

Jk/2 J\q—k\ Z 

+ % gh (q jP , k)P gf (p)P hb (k) + j fgh (p : k, q)P gb (k)P he (q)} + (q <-> p)} . (B.6) 

Comparing (IB.ip to flB.6p . we see that the latter integral contains only running functions 
of one variable, the P's. 

The evolution of the PS can then be obtained by solving the system of eqs. (IB. 41) 
and (1B.5H . By inspection of the explicit expression for the vertex, eq. ([HD, and of the 
symmetry properties of the integral (1B.6|) . one can conclude that of the 64 components 
of I a cd,bef{k) only 12 take part to the system. Indeed, the 12 independent components 
to follow are identified by the direct product of the two triplets (acd) = (112), (222) and 
the six triplets (fee/) = (fell), (fel2), (622), (6 = 1, 2). Together with the 3 equations for 
the independent components of the PS (namely, Pn, Pu, and P22), we have a total of 
15 equations. 

It is also convenient to put the integrals such as eq. (IB.6P in the following form, 

poo rq 

/ dqq dpp \F(k,q,p) + (q <-* p)\ 

Jk/2 J\a-k\ 



Ik/2 -]\q-k\ 

roo rk/s/2 X 2 — V 2 

I dx dy ■ 

Jk/V2 JO 



T?n x + y x-y 



,(B-7) 



lk/V2 JO 2 

which has the virtue of having both integration intervals independent on the integration 
variables. 
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